rm(list=ls())
graphics.off()

library(readxl)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(ggpubr)
library(rstatix)
library(nlme)
library(matrixcalc)
library(Biodem)
library(stats)
library(multcomp)
library(hrbrthemes)
library(viridis)
library(gapminder)
library(devtools)
library(stringr)
library(patchwork)
library(gridExtra)
library(cowplot)
library(extrafont)
library(gtable)
library(ggtext)
library(effsize)
library(data.table)
library(psych)
library(psychotools)
library(psychTools)
library(BlandAltmanLeh)

######################### START ###########################################

setwd("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3")

data <- read.csv("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3/IMU_Mean.csv")


#Tables
AC <- bind_cols(data[,1],data[,2:7])
AT <- bind_cols(data[,1],data[,8:13])


colnames(AC) <- c("ID","aclow","bcchall","cchigh","dmlow","emchall","fmhigh")
colnames(AT) <- c("ID","aclow","bcchall","cchigh","dmlow","emchall","fmhigh")

###AC 
AC_long <- pivot_longer(AC, cols=c("aclow","bcchall","cchigh","dmlow","emchall","fmhigh"))
div <- rep(c("Mental","Motor"), each = 3)
AC_long$div <- rep(div,nrow(AC))

AC_long = split(AC_long, AC_long$div)

AC_coglong = AC_long$Mental[,1:3]

AC_motlong =  AC_long$Motor[,1:3]

###

plot <- ggplot(AC_long, aes(name, value, ID, fill = name)) +
  geom_boxplot() +
  labs(x = 'Intensity Level',
       title = 'AC',
       y = "AC") +
  theme(plot.title = element_text(hjust = 0.5, size = 20),
        legend.title=element_blank(),
        panel.background = element_blank(),
        panel.grid.major.y = element_line(colour = "grey", size = 0.5),
        axis.title = element_text(size = 13),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 13)) +
  scale_x_discrete(label = element_blank()) +
  scale_y_continuous(limits = c(-0.015,0.223)) +
  guides(fill="none") +
  scale_fill_manual(values=c("steelblue4", "lightblue1", "cadetblue3", "steelblue4", "lightblue1", "cadetblue3")) +
  annotate(geom = "text", x = 0.83, y = -0.009, label = "Very Low", hjust = 0, vjust = 0, size = 3, color = "gray50") +  
  annotate(geom = "text", x = 0.85, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 1.85, y = -0.009, label = "Challenging", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 1.85, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 2.85, y = -0.009, label = "Very High", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 2.87, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  facet_grid(. ~ div, scales = "free", space = "free")


#ggsave("HRVplot.png", plot, width = 10.5, height = 6.92)

###Statistics

###Test Normality
Norm.AC.CogEasy = shapiro.test(AC$aclow)
Norm.AC.CogMiddle = shapiro.test(AC$bcchall)
Norm.AC.CogDifficult = shapiro.test(AC$cchigh)

Norm.AC.MotEasy = shapiro.test(AC$dmlow)
Norm.AC.MotMiddle = shapiro.test(AC$emchall)
Norm.AC.MotDifficult = shapiro.test(AC$fmhigh)

###Perform Test and Post-Hoc
### Cognitive 
if (Norm.AC.CogEasy$p.value > 0.05 & Norm.AC.CogMiddle$p.value > 0.05 & Norm.AC.CogDifficult$p.value > 0.05){
  print(anova_test(data=AC_coglong, dv = value, wid=ID, within = name))
  pairwise.t.test(AC_coglong$value, AC_coglong$name, p.adjust.method = 'bonferroni', paired = TRUE)
} else { 
  print(friedman.test(data=AC_coglong, value ~ name|ID))
  pairwise.wilcox.test(AC_coglong$value, AC_coglong$name, p.adjust.method = 'bonferroni', paired = TRUE)
}

###Motor
### Cognitive 
if (Norm.AC.MotEasy$p.value > 0.05 & Norm.AC.MotMiddle$p.value > 0.05 & Norm.AC.MotDifficult$p.value > 0.05){
  print(anova_test(data=AC_motlong, dv = value, wid=ID, within = name))
  pairwise.t.test(AC_motlong$value, AC_motlong$name, p.adjust.method = 'bonferroni', paired = TRUE)
} else { 
  print(friedman.test(data=AC_motlong, value ~ name|ID))
  pairwise.wilcox.test(AC_motlong$value, AC_motlong$name, p.adjust.method = 'bonferroni', paired = TRUE)
}


#ICC analysis

data_AC_Test <- read.csv("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3/IMU_Day1.csv")
data_AC_ReTest <- read.csv("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3/IMU_Day2.csv")


## AC

## Cognitive

### Easy
data_AC_Cog_easy = data.table(data_AC_Test$AC_CognitionEasy, data_AC_ReTest$AC_CognitionEasy)
#data_AC_easy_RRn = data_AC_easy_RRn[-3,] #delete unwanted ID (1, 12, 14?)

ICC_AC_Cog_easy =ICC(data_AC_Cog_easy, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_AC_Cog_easy<- bland.altman.plot(data_AC_Test$AC_CognitionEasy,data_AC_ReTest$AC_CognitionEasy, main="Test-Retest AC Easy", xlab="Means", ylab="differences")

print(ICC_AC_Cog_easy)

###Middle
data_AC_Cog_middle = data.table(data_AC_Test$AC_CognitionOptimal, data_AC_ReTest$AC_CognitionOptimal)
#data_AC_Cog_middle = data_AC_Cog_middle[-3,] #delete unwanted ID (1, 12, 14?)

ICC_AC_Cog_middle =ICC(data_AC_Cog_middle, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_AC_Cog_middle <-bland.altman.plot(data_AC_Test$AC_CognitionOptimal,data_AC_ReTest$AC_CognitionOptimal, main="Test-Retest AC Optimal", xlab="Means", ylab="differences")

print(ICC_AC_Cog_middle)


###Difficult
data_AC_Cog_difficult = data.table(data_AC_Test$AC_CognitionDifficult, data_AC_ReTest$AC_CognitionDifficult)
#data_AC_Cog_difficult = data_AC_Cog_difficult[-3,] #delete unwanted ID (1, 12, 14?)

ICC_AC_Cog_difficult =ICC(data_AC_Cog_difficult, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_AC_Cog_difficult<- bland.altman.plot(data_AC_Test$AC_CognitionDifficult,data_AC_ReTest$AC_CognitionDifficult, main="Test-Retest AC Difficult", xlab="Means", ylab="differences")

print(ICC_AC_Cog_difficult)

## Motor

### Easy
data_AC_Mot_easy = data.table(data_AC_Test$AC_MotorEasy, data_AC_ReTest$AC_MotorEasy)
#data_AC_easy_RRn = data_AC_easy_RRn[-3,] #delete unwanted ID (1, 12, 14?)

ICC_AC_Mot_easy =ICC(data_AC_Mot_easy, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_AC_Mot_easy<- bland.altman.plot(data_AC_Test$AC_MotorEasy,data_AC_ReTest$AC_MotorEasy, main="Test-Retest AC Easy", xlab="Means", ylab="differences")

print(ICC_AC_Mot_easy)

###Middle
data_AC_Mot_middle = data.table(data_AC_Test$AC_MotorOptimal, data_AC_ReTest$AC_MotorOptimal)
#data_AC_Mot_middle = data_AC_Mot_middle[-3,] #delete unwanted ID (1, 12, 14?)

ICC_AC_Mot_middle =ICC(data_AC_Mot_middle, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_AC_Mot_middle <-bland.altman.plot(data_AC_Test$AC_MotorOptimal,data_AC_ReTest$AC_MotorOptimal, main="Test-Retest AC Optimal", xlab="Means", ylab="differences")

print(ICC_AC_Mot_middle)


###Difficult
data_AC_Mot_difficult = data.table(data_AC_Test$AC_MotorDifficult, data_AC_ReTest$AC_MotorDifficult)
#data_AC_Mot_difficult = data_AC_Mot_difficult[-3,] #delete unwanted ID (1, 12, 14?)

ICC_AC_Mot_difficult =ICC(data_AC_Mot_difficult, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_AC_Mot_difficult<- bland.altman.plot(data_AC_Test$AC_MotorDifficult,data_AC_ReTest$AC_MotorDifficult, main="Test-Retest AC Difficult", xlab="Means", ylab="differences")

print(ICC_AC_Mot_difficult)



###################################################################################################################################################


###AT 
AT_long <- pivot_longer(AT, cols=c("aclow","bcchall","cchigh","dmlow","emchall","fmhigh"))
div <- rep(c("Mental","Motor"), each = 3)
AT_long$div <- rep(div,nrow(AT)) 

AT_long = split(AT_long, AT_long$div)

AT_coglong = AT_long$Mental[,1:3]

AT_motlong =  AT_long$Motor[,1:3]

###

plot <- ggplot(AT_long, aes(name, value, ID, fill = name)) +
  geom_boxplot() +
  labs(x = 'Intensity Level',
       title = 'AT',
       y = "AT") +
  theme(plot.title = element_text(hjust = 0.5, size = 20),
        legend.title=element_blank(),
        panel.background = element_blank(),
        panel.grid.major.y = element_line(colour = "grey", size = 0.5),
        axis.title = element_text(size = 13),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        strip.text.x = element_text(size = 13)) +
  scale_x_discrete(label = element_blank()) +
  scale_y_continuous(limits = c(-0.015,0.223)) +
  guides(fill="none") +
  scale_fill_manual(values=c("steelblue4", "lightblue1", "cadetblue3", "steelblue4", "lightblue1", "cadetblue3")) +
  annotate(geom = "text", x = 0.83, y = -0.009, label = "Very Low", hjust = 0, vjust = 0, size = 3, color = "gray50") +  
  annotate(geom = "text", x = 0.85, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 1.85, y = -0.009, label = "Challenging", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 1.85, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 2.85, y = -0.009, label = "Very High", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  annotate(geom = "text", x = 2.87, y = -0.015, label = "Intensity", hjust = 0, vjust = 0, size = 3, color = "gray50") +
  facet_grid(. ~ div, scales = "free", space = "free")


#ggsave("HRVplot.png", plot, width = 10.5, height = 6.92)

###Statistics

###Test Normality
Norm.AT.CogEasy = shapiro.test(AT$aclow)
Norm.AT.CogMiddle = shapiro.test(AT$bcchall)
Norm.AT.CogDifficult = shapiro.test(AT$cchigh)

Norm.AT.MotEasy = shapiro.test(AT$dmlow)
Norm.AT.MotMiddle = shapiro.test(AT$emchall)
Norm.AT.MotDifficult = shapiro.test(AT$fmhigh)

###Perform Test and Post-Hoc
### Cognitive 
if (Norm.AT.CogEasy$p.value > 0.05 & Norm.AT.CogMiddle$p.value > 0.05 & Norm.AT.CogDifficult$p.value > 0.05){
  print(anova_test(data=AT_coglong, dv = value, wid=ID, within = name))
  pairwise.t.test(AT_coglong$value, AT_coglong$name, p.adjust.method = 'bonferroni', paired = TRUE)
} else { 
  print(friedman.test(data=AT_coglong, value ~ name|ID))
  pairwise.wilcox.test(AT_coglong$value, AT_coglong$name, p.adjust.method = 'bonferroni', paired = TRUE)
}

###Motor
### Cognitive 
if (Norm.AT.MotEasy$p.value > 0.05 & Norm.AT.MotMiddle$p.value > 0.05 & Norm.AT.MotDifficult$p.value > 0.05){
  print(anova_test(data=AT_motlong, dv = value, wid=ID, within = name))
  pairwise.t.test(AT_motlong$value, AT_motlong$name, p.adjust.method = 'bonferroni', paired = TRUE)
} else { 
  print(friedman.test(data=AT_motlong, value ~ name|ID))
  pairwise.wilcox.test(AT_motlong$value, AT_motlong$name, p.adjust.method = 'bonferroni', paired = TRUE)
}


#ICC analysis

data_AT_Test <- read.csv("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3/IMU_Day1.csv")
data_AT_ReTest <- read.csv("//nausers01/User/GOIKOG/Desktop/MatLab/Output/IntensityMeasures3/IMU_Day2.csv")


## AT

## Cognitive

### Easy
data_AT_Cog_easy = data.table(data_AT_Test$AT_CognitionEasy, data_AT_ReTest$AT_CognitionEasy)
#data_AT_easy_RRn = data_AT_easy_RRn[-3,] #delete unwanted ID (1, 12, 14?)

ICC_AT_Cog_easy =ICC(data_AT_Cog_easy, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_AT_Cog_easy<- bland.altman.plot(data_AT_Test$AT_CognitionEasy,data_AT_ReTest$AT_CognitionEasy, main="Test-Retest AT Easy", xlab="Means", ylab="differences")

print(ICC_AT_Cog_easy)

###Middle
data_AT_Cog_middle = data.table(data_AT_Test$AT_CognitionOptimal, data_AT_ReTest$AT_CognitionOptimal)
#data_AT_Cog_middle = data_AT_Cog_middle[-3,] #delete unwanted ID (1, 12, 14?)

ICC_AT_Cog_middle =ICC(data_AT_Cog_middle, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_AT_Cog_middle <-bland.altman.plot(data_AT_Test$AT_CognitionOptimal,data_AT_ReTest$AT_CognitionOptimal, main="Test-Retest AT Optimal", xlab="Means", ylab="differences")

print(ICC_AT_Cog_middle)


###Difficult
data_AT_Cog_difficult = data.table(data_AT_Test$AT_CognitionDifficult, data_AT_ReTest$AT_CognitionDifficult)
#data_AT_Cog_difficult = data_AT_Cog_difficult[-3,] #delete unwanted ID (1, 12, 14?)

ICC_AT_Cog_difficult =ICC(data_AT_Cog_difficult, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_AT_Cog_difficult<- bland.altman.plot(data_AT_Test$AT_CognitionDifficult,data_AT_ReTest$AT_CognitionDifficult, main="Test-Retest AT Difficult", xlab="Means", ylab="differences")

print(ICC_AT_Cog_difficult)

## Motor

### Easy
data_AT_Mot_easy = data.table(data_AT_Test$AT_MotorEasy, data_AT_ReTest$AT_MotorEasy)
#data_AT_easy_RRn = data_AT_easy_RRn[-3,] #delete unwanted ID (1, 12, 14?)

ICC_AT_Mot_easy =ICC(data_AT_Mot_easy, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_AT_Mot_easy<- bland.altman.plot(data_AT_Test$AT_MotorEasy,data_AT_ReTest$AT_MotorEasy, main="Test-Retest AT Easy", xlab="Means", ylab="differences")

print(ICC_AT_Mot_easy)

###Middle
data_AT_Mot_middle = data.table(data_AT_Test$AT_MotorOptimal, data_AT_ReTest$AT_MotorOptimal)
#data_AT_Mot_middle = data_AT_Mot_middle[-3,] #delete unwanted ID (1, 12, 14?)

ICC_AT_Mot_middle =ICC(data_AT_Mot_middle, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_AT_Mot_middle <-bland.altman.plot(data_AT_Test$AT_MotorOptimal,data_AT_ReTest$AT_MotorOptimal, main="Test-Retest AT Optimal", xlab="Means", ylab="differences")

print(ICC_AT_Mot_middle)


###Difficult
data_AT_Mot_difficult = data.table(data_AT_Test$AT_MotorDifficult, data_AT_ReTest$AT_MotorDifficult)
#data_AT_Mot_difficult = data_AT_Mot_difficult[-3,] #delete unwanted ID (1, 12, 14?)

ICC_AT_Mot_difficult =ICC(data_AT_Mot_difficult, missing=TRUE, alpha = 0.05, lmer = TRUE, check.keys = FALSE) #last 2 variables not needed, they are so by default

BAP_AT_Mot_difficult<- bland.altman.plot(data_AT_Test$AT_MotorDifficult,data_AT_ReTest$AT_MotorDifficult, main="Test-Retest AT Difficult", xlab="Means", ylab="differences")

print(ICC_AT_Mot_difficult)